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Abstract. We present global 2D and 3D simulations of self-gravitating 

\^ ' magnetized tori. We used the 2D calculations to demonstrate that the 

. properties of the MRI are not affected by the presence of self-gravity: 

^ ' MHD turbulence and enhanced angular momentum transport follow the 

^ . linear growth of the instability. In 3D, we have studied the interaction 

<^ ' between an m = 2 gravitational instability and MHD turbulence. We 

. found its strength to be significantly decreased by the presence of the 

I latter, showing that both instabilities strongly interact in their non-linear 

^ ■ phases. We discuss the consequences of these results. 

o ■ 
o : 

+l3 , 1. Introduction 

. 

^ I In the early phases of star formation, the forming accretion disks are likely to 

•'-j ■ be very massive because of a strong infall from the parent molecular cloud. 

rS I These massive disks are subject to the development of gravitational instabilities 

■ which transport angular momentum outward (Laughlin, Korchagin, & Adams 

1998). In addition, when sufficiently ionized, the disks are also unstable to the 
MRI (Balbus & Hawley 1991, 1998). The simultaneous development of these 
instabilities in disks may significantly affect their evolution. We have undertaken 
a study of self-gravitating magnetized tori by means of numerical simulations. 

In section 2, we describe the numerical methods we have used. In section 
3, we present the results of the 2D simulations, focusing on the properties of 
the MRI in a self-gravitating environment. In section 4, we review preliminary 
results obtained in 3D calculations, and we discuss the implications of our results 
in section 5. 
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2. Numerical methods 

We used the code Zeus-2D (Stone Sz Norman 1992a, 1992b) to perform the 
2D calculations and the code GLOBAL (Hawlcy &: Stone 1995) for the 3D 
simulations. Both codes solve the MHD equations using time-explicit Eulcrian 
finite differencing and the Constrained Transport method to evolve the magnetic 
field. The Poisson solver used to calculate the gravitational potential $g involves 
two steps: we first calculate on the boundary, using the Legendre functions 
well-suited to the cylindrical geometry (Cohl 1999), and we apply the Successive 
Over-Relaxation method to update everywhere on the grid (Hirsch 1988). 

3. 2D simulations 

We have built an initial equilibrium structure using the Self— Consistent Field 
(SCF) method (Hachisu 1986). The inner and outer radii of the torus arc = 
0.3 and Rout = 1 respectively. We set the angular velocity O oc r~^'^^, where 
r is the cylindrical radius. A central mass having half the mass of the disk is 
present. The equation of state is adiabatic. 

In the following, we compare the evolution of this torus to its zero mass 
counterpart, the gravitational potential of the latter being only that of a central 
point mass. At the beginning of the simulation, a weak poloidal magnetic field 
is added to the equilibrium structure, with the toroidal component of the vector 
potential being: 



where p is the mass density. The components of the magnetic field are then scaled 
such that the volume averaged ratio of gas to magnetic pressure (hereafter called 
(/?)) equals 1500 for the self-gravitating model and 200 for the zero mass one. 
The resolution in both model is {Nr,Nz) = (256, 256) in the radial and vertical 
directions, respectively. 

We found the MRI grows in both models, developing approximately the 
same Maxwell stress (defined by T^"-^ = —BrB^/An). In the self-gravitating 
case, the evolution is very similar to what was found in previous non self- 
gravitating calculations (e.g. Hawley [2000)]. Namely, the initial linear growth 
of the instability is followed by a turbulent phase during which angular momen- 
tum is transported outward. Turbulence then gradually decays because of the 
anti-dynamo theorem. 

We show in figure 1 the logarithm of the density field in the (r — z) plane 
during the turbulent phase for both models. In the self-gravitating case {left 
panel), the initial torus has developed a two-component structure, composed of 
an inner thin disk fed by an outer thick, massive torus. Figure 2 shows the 
angular momentum I radial profile in the equatorial plane during this phase 
{solid line). As shown by the dashed line, the inner disk is in Keplerian rotation 
around the central mass. The dotted line is a fit of the outer part of the disk 
with a power law dependence I oc r°-^, very close to the Mestel profile I oc r. 

Coming back to figure 1, the comparison between both models is striking: 
although the Maxwell stress is similar in the two cases, the self-gravitating torus 
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Figure 1. Comparison between the logarithm of the density in the 
(r-z) plane for the self-gravitating torus {left panet) and the zero mass 
torus (right panel). The former develops a two-component structure 
composed of an inner Keplerian disks and an outer, more massive, 
thick disk. The later is disrupted by the MRI, and a standard thin 
disk structure builds up, with a constant H/R ratio. 



presents a much more featureless structure than its zero mass counterpart. This 
result is due to the large scale coherence of the self-gravitating potential, which 
smoothes MRI fluctuations. We shall see that this effect persists in 3D. 



4. 3D simulations 



In this section, we present the results of our 3D simulations. As above, we have 
used the SCF method to compute the disk model, but with different parame- 
ters. Here the central mass is twice that of the disk, and the specific angular 
momentum / in the disk is chosen such that 



l(xn?{r), (2) 

where m(r) is the mass inside the radius r divided by the total disk mass. 
The density contours of this disk model are shown in figure 3 (left panel). We 
introduce the Toomre parameter Q defined by: 

where Cg is the sound speed, k the epicyclic frequency and S the surface density 
of the disk. Linear studies have shown that a disk becomes gravitationally 
unstable when Q is of the order of 1. The radial profile of Q in our disk model 
is shown on the right panel of figure 3. We have Q ^ 1 over a large range of 
radii. We therefore expect the disk to be gravitationally unstable. 

4.1. Hydrodynamical simulations 

We have performed hydrodynamical simulations to check that our code behaves 
as expected in this regime. In the simulations presented here, the resolution is 
(Nr, Nfp, Nz) = (128,64,64) in cylindrical coordinates. The azimuthal domain 
is restricted to ^ G [0, tt], which prevents the growth of modes with odd values 
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Figure 2. Angular momentum profile in the equatorial plane of the 
torus, averaged between orbits 5.05 and 7. The dashed line shows the 
Keplerian profile resulting from the central point mass. The dotted 
line has a power law dependence with radius: I cx r^-^. 



of the azimuthal number m. Time is measured in units of the orbital period at 
the initial outer edge of the disk. 

The time evolution of the Fourier component of the density in the equatorial 
plane is shown in figure 4 for the modes m = 2,4,6,8 (from top to bottom). 
The m = 2 mode grows at the beginning of the simulation and saturates after 4 
orbits. The development of this 2 arms spiral structure can be seen on figure 5, 
which shows the logarithm of the density in the equatorial plane after 4.27 orbits. 

During the simulation, matter is driven toward the center of the disk by the 
gravitational torque associated with the spiral arms. As the surface density in 
the disk decreases, the Toomre Q parameter increases. As a result, the instability 
dies away after 8 orbits. 

These results are in agreement with earlier simulations of similar disks. 

4.2. 3D MHD simulations in an axisymmetric potential 

Next, we have performed a 3D MHD calculation of the evolution of the disk in an 
axisymmetric gravitational potential, i.e. a potential for which wc have retained 
only the m = component. This prevents the growth of gravitational insta- 
bilities and enables us to check the behavior of our code in the MHD regime. 
The resolution used for this run is {Nr, N^, N^) = (128,32,128) and the az- 
imuthal domain is [0,7r/2]. A toroidal magnetic field with = 8 is added to 
the equilibrium model. 

The left panel of figure 6 shows the time evolution of the volume averaged 
Maxwell and Reynolds stress tensors. The Maxwell stress rises during the early 
part of the evolution (corresponding to the linear development of the instability) 
before the MRI breaks down into turbulence after 4 orbits. The Reynolds stress 
is smaller during the whole simulation. The right panel of figure 6 shows the 
radial profile of the a parameter (defined as the ratio of the total stress to the 
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Figure 3. Density contours for the disk model used in the 
3D simulations (left panel). The contours shown arc p = 
lO"'^, 10"^, 0.001, 0.01, 0.1, 0.5, 0.7 and 0.9. The right panel shows the 
radial profile of the Toomre Q parameter. 
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Figure 4. Time evolution of the Fourier components of the density in 
the equatorial plane for the modes m = 2, 4, 6, 8 (from top to bottom). 
The m = 2 mode is unstable. 

thermal pressure) at the end of the simulation. It has a typical value of a few 
times 10~^, similar to values seen in 3D simulations of zero mass disks starting 
with similar magnetic field configurations (Steinacker Sz Papaloizou 2002). 

4.3. Full 3D simulations 

In this section, we present the results of the simulation in which both the MRI 
and gravitational instabilities can develop. We include the Fourier components 
of the gravitational potential up to m = 8. The toroidal magnetic field that we 
add to the equilibrium structure is the same as that described in the previous 
section. For this calculation, the resolution is {Nr,N^,Nz) = (128,64,128) and 
the azimuthal domain is [0,7r]. 

As expected, we observe the simultaneous appearance of MHD turbulence 
along with the m = 2 spiral arm reminiscent of the hydrodynamical calculation. 
To analyze the angular momentum transport properties, we compare the time 
evolution of the volume averaged gravitational and Maxwell stress tensors ob- 
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Figure 5. Logarithm of the density in the equatorial plane after 4.27 
orbits. 

tained in this simulation with that obtained in the runs described above. Note 
that the gravitational stress tensor is defined by: 

~ AttgJ r dr d<i> ^ ' 

The left panel of figure 7 shows ™^ as a function of time for this simulation 
{dotted line) and during the hydrodynamical run {solid line). We can see that the 
presence of MHD turbulence decreases the strength of the gravitational stress by 
a factor of about 2. This result shows that MHD turbulence and the gravitational 
instability strongly interact with each other. The time evolution of the Maxwell 
stress is shown on the right panel of figure 7 for the full simulation {dotted line) 
and for the MHD run in an axisymmetric potential {dashed line). We see that 
the Maxwell stress increases when the gravitational instability develops. This 
appears to be due to the accumulation of the magnetic field along the spiral arms. 
When the gravitational instability disappears after about 7 orbits (because of a 
decrease of the mass density), the Maxwell stress goes back to its "turbulent" 
value. 

To test the sensitivity of our results to the initial field topology, we have 
run the same calculation with an initial poloidal magnetic field. We found the 
same qualitative results. 

5. Conclusions and Perspectives 

We have presented here the first global simulations of massive, magnetized disks. 

Using 2D axisymmetric numerical simulations, we have shown that the 
properties of the MRI are similar in self-gravitating and zero mass disks. We 
observe that these disks quickly develop a dual structure composed of an inner 
thin disk in Keplerian rotation around the central mass, and a thicker outer 
torus whose rotation profile is close to a Mestel profile. 

We have then used 3D simulations to study the angular momentum trans- 
port properties in disks when both MHD turbulence and gravitational instabil- 
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Figure 6. Left panel: time evolution of the Maxwell {solid line) and 
Reynolds (dashed line) stress tensor in the MHD run with an axisym- 
metric gravitational potential. Right panel: radial profile of the a 
parameter. 



itics arc present. Wc have found that the gravitational instability is affected by 
the presence of the turbulence: the gravitational stress tensor is decreased by 
roughly a factor of 2 when compared with hydrodynamical simulations. This re- 
sults in a smaller mass accretion rate toward the central object. Self-gravitating 
disks may therefore have a longer lifetime than previously thought. 

Note that the simulations presented here use an adiabatic equation of state. 
In this case, all the energy generated in shocks and compression is locally con- 
verted into heat. This prevents the formation of bound objects by gravitational 
collapse. The opposite case would correspond to the use of a locally isothermal 
equation of state, for which all the energy generated is immediately radiated 
away. Several authors (Mayer et al. 2002, Rice et al. 2003, Boss 1997 k 1998) 
have indeed reported gravitational collapse in isothermal calculations of disks, 
although this issue is still under debate (Pickett et al. 2003). We beheve that 
the presence of MHD turbulence affects the thermal balance in the disk and 
therefore needs to be included. We are currently performing calculations of 
isothermal massive and magnetized disks. 
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